Simulation of granular jet: Is granular flow really a "perfect fluid? 
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We perform a three-dimensional simulation of granular jet for both frictional and frictionless 
grains. Through our simulation we confirm that the system after an impact of the jet has a finite 
viscosity for both frictional and frictionless grains. This result is contrast to the suggestion in an 
experiment of granular jet [X. Cheng et at, Phys. Rev. Lett. 99, 188001 (2007) ], in which there 
exists an analogy between granular jets and quark-gluon plasma. 
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PACS numbers: 45.70.Mg, 83.10.Rs,05.20.Dd 

Introduction: The impact processes play important 
roles in various fields such as nuclear reactions [H-Q, 
atomic collisions d, hydrodynamics and granular 
physics [j]HI3. 

Recent experimental and numerical stud- 
ies revealed interesting aspects of impact processes of a 
granular flow. At low volume fractions, impact of a gran- 
ular flow onto a wall produces a shock, which quantita- 
tively a gree s with the Mach cone produced by supersonic 
gas fk>w [lll4l3j . The impact dynamics of granular parti- 
cles is important not only for industrial app lications, e.g. 
ink-jet printing and blast cleaning [13, Il8| . but also for 
geophysical problems such as the formation of craters [13- 

Recently, an experimental paper on dense granular 
jets[l(| has reported that fluid state after the impact is 
similar to that for quark-gluon plasma (QGP) achieved 
in heavy ion hadron colliders, where QGP behaves as a 
fluid with very small viscosity QHH. Quite recently, El- 
lowitz et al. [19|] demonstrated that the solution of inviscid 
Euler equation is almost identical to that obtained from 
a molecular dynamics simulation for inelastic hard core 
particles, at least, for two dimensional frictionless grains. 
These results are counter intuitive because the dense 
granular fluid has a large viscosity in usual setup [20]. 

The purpose of this Letter is to clarify whether the 
granular fluid after the jet impact on a fixed wall ac- 
tually behaves as a perfect fluid. So far all numerical 
studiesfl9l.l2ll.l22l] are two-dimensional ones for this prob- 
lem, we perform a three dimensional molecular dynamics 
simulation for soft core particles to compare our results 
with experimental results. 

Model: We adopt the discrete element method (DEM) 
for mono-disperse soft core particles of the diameter 
[24|. The reason why we adopt a soft core model 
is the following. The DEM can be used even for dense 
systems above the jamming point, while the simulation 
in terms of the event driven algorithm cannot be free 
from inelastic collapse [25jj and cannot reach the jamming 
point. DEM has also an advantage when we consider the 
effect of friction and rotation of grains. Indeed, the event- 
driven algorithm for frictional grains is complicated, but 
DEM for frictional grains is simple. When the particle i 
at the position and the particle j at rv,- are in contact, 
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FIG. 1: (Color online) Snapshot of a three-dimensional sim- 
ulation. 



nj = |r, — Tj\ and gy = Vj — Vj with the velocity v» of 
the particle i. The tangential contact force is given by 
Fh = rmn{|Fy|,/ii ? i "}sign(-F' i *.), where fi is the friction 
constant, sign(x) = 1 for x > and sign(a;) = — 1 for oth- 
erwise, Ffj = ktSjj — rjtSjj with the tangential overlap Sjj 
between i and j particles and the tangential component of 
relative velocity Sfj between i th and j th particles. Here, 
we adopt parameters k t — 0.2k ni rj t — 0.5ry„, /j = 0.2, 
k n = 4.98 x I0 2 mul/d 2 and r/„ = 2.88uo/d, where uo 
is the incident velocity and m is the particle mass. Ex- 
perimentally, it is known that the friction constant of 
spheres is fx = 0.175 ± 0.1 [2^] . This set of parameters 
implies that the restitution coefficient for normal impact 
is e = 0.75 and duration time ist c = 0.10d/uo. We adopt 
the second-order Adams-Bashforth method for the time 
integration with the time interval At = 0.02i c through 
our simulation. We fix the jet diameter D- ]C t/d = 10.0 
and the target diameter D tar /d = 22.0. We generate an 
initial configuration as follows. We prepare fee crystals 
and remove particles randomly to reach the desired den- 
sity. The initial granular temperature, which represents 
the fluctuation of particles motion, is zero. We control 
the initial volume fraction 4>o/4>i cc = (j)Q before the im- 
pact as 0.30 < cf>o < 0.90 with volume fraction for a fee 
crystal </)f cc ~ 0.74. We use 20,000 particles for our sim- 
ulations. The wall consists of one-layer particles, which 
are connected each other and with their own initial po- 
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sitions via the spring and the dashpot. Here, the spring 
constant and the dashpot constant for the wall particles 
are given by k p — IO-Otouq/cP and r] p — 5.0ry„, respec- 
tively. It is known that the collective motion of particles 
near the wall is almost frozen [19J, while the granular tem- 
perature near the wall is higher than that in the other 
regions. 

Figure 1 is a snapshot of our simulation on the im- 
pact of granular jet. The scattering angle exhibits the 
crossover cone-like structure to sheet-like one, depending 
on D tar / Dj et , which is almost independent of the resti- 
tution coefficient [13]. 

We evaluate physical quantities near the wall whose 
height is z — Az = 5.0d from the wall z = in our 
simulations. We divide the cylindrical calculation re- 
gion into the radial direction r = 0, Ar, 2Ar, • • • , 5Ar, 
with Ar = i?tar/5 and the target radius i?t a r- Phys- 
ical quantities are estimated in k th mesh region with 
kAr <r <(k + l)Ar (k = 0, 1, • • • , 5) and < z < Az. 

Stress tensor. Let us evaluate the stress tensor near 
the wall as in Ref. 28]. The stress tensor (T a ^(r) at r 
consists of the kinetic part cr^ (r) and the contact stress 
<J c a/3 {r) as 



o- a /3(r) = <J k ap {r) +a c ap (r), 



(1) 



where their microscopic definitions are respectively given 
by 
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where i and j are indices of particles, a, j3 = r,8, z de- 
notes cylindrical coordinates and denotes the sum- 
mation over the particles located at r. Here, z axis is 
parallel to the incident jet axis, and V is the volume of 
each mesh at r and Ui a = v l a — v a with the mean veloc- 
ity v a . To calculate the stress tensor in the cylindrical 
coordinates, we firstly calculate <j a '$' in the Cartesian 
coordinate, a',/3' = x,y,z, whose original point is the 
same as the cylindrical one, and transform it into that 
for the cylindrical one. 

Here we show the profile of the stress tensor for the fric- 
tional case (Fig. [5] for (f> — 0.90). We average the data 
over ten different initial configurations. From Fig. [3J 
it is apparent that off-diagonal components of the stress 
tensor a rz and <r zr are much smaller than diagonal com- 
ponents, where the ratio of the off-diagonal to the di- 
agonal element is estimated as \<J rz ja zz \ ~ 1.7 x 10~ 2 
at r/Rtm = 0.1. This result supports that the solution 
of Euler equation well reproduce the granular flow af- 
ter the impact 19]. We also found that there exists a 
large normal stress difference, which cannot be observed 
even for our two dimensional case. We obtain the ratio 
Wrz/o- zz \ ~ 3.0 x 10 -2 at r/i?tar = 0.1 for the frictionless 
case, where off-diagonal components are much smaller 
than diagonal ones as in the case of the frictional case. 
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FIG. 2: (Color online) The profile of the stress tensor as the 
functions of distance from the jet axis r with R ta , r for fric- 
tional grains with (j>o = 0.90. The off-diagonal components of 
the stress tensor oy z and o zr are much smaller than diagonal 
components as \a rz /a zz \ ~ 1.7 x 10 -2 at r/i? t ar = 0.1. 



Pressure: Let us look at the result of the pressure here 
(Fig. [3]). It is known that the granular sheared flow 
such as a chute flow and a plane shear flow can be ap- 
proximately described by granular hydrodynamics with 
transport coefficient derived from the kinetic theory, i.e. 
the Enskog equation [29j - [3^ ] . We adopt the transport co- 
efficients derived by Garzo and Duftyj!o[ with the pres- 
sure P = <7 aa /3, the density n, the volume fraction 



4>, and the granular temperature T g (r) = 
The pressure is conventionally given by 

— = 1 + 20(1 + e) X , 
nl„ 



ma 
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where <j)f — 0.49 and (f> c — 0.64[35|. For the frictional 
case, in general, five equations for rotational degree of 
freedom are necessary, in addition to those for transla- 
tional one. However, it is known that ten equations for 
frictional grains can be reduced to five equations by intro- 
ducing effective coefficient of restitution e if the friction 
constant /i is small [3214341 ] . According to this simplifica- 
tion we use the effective restitution coefficient e = 0.616 
for e = 0.75 and /i = 0.2, for frictional case in the follow- 
ing. 

Let us compare the theoretical curve with numerical 
results for several <fio (Fig. [3]). The black solid line in 
Fig. [3] and that in the inset denote the theoretical curve 
for the frictionless case and the frictional case, respec- 
tively. Surprisingly, the expression for the pressure in 
Eq. ^ considerably reproduces the numerical result for 
<j> < 0.5 inspite of the existence of the normal stress dif- 
ference, while Eq.Q for 0.5 < <fi < 0.6 may have signifi- 
cant deviation from the theoretical line. The deviation, 
which may result from the singularity at the source point 
around r ~ 0, emerges only at r / R tav = 0.1. 
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FIG. 3: (Color online) The comparison between the theoreti- 
cal pressure in Eq. (4) and the observed pressure for several 
4>q, where the vertical axis is P divided by the density n and 
the granular temperature T g for the frictionless case. The in- 
set denotes comparison between those for the frictional case. 
Black solid lines in each figures denote Eq. (4) for the fric- 
tionless and the frictional case, respectively. 
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FIG. 4: (Color online) The diagonal components of the 
stress tensor for the frictional case with several <j>o di- 
vided by nT a , where T a is the temperature for a direc- 
tion T a = ^2 i mUi cc /N. Red empty points, blue filled 
points, purple half-filled points and the solid black line de- 
note P z /nT z , P r /nT r , Pg/nTg and Eq. (4), respectively. 



Although there exists large normal stress differences, 
numerical results can be reproduced from the kinetic the- 
ory if we introduce the anisotropic temperature. Indeed, 
equations of state for each coordinates satisfies, P a = 
nT a {l+20(l+e)x} for a = r,9, z and P r — a rr , Pg = ogg 
and P z — cr zz . By summing up P a = nT a {l + 20(1 + e)x} 
over a, we obtain Eq. (4). From Fig. 4, we verify 
that P z /nT z and P r /nT r are on the theoretical curve 
for isotropic systems, but Pg/nTg has a systematically 
larger value from the isotropic one. Although our sugges- 
tion that the anisotropy only appears through the normal 
stress is not perfect, the result gives a reasonable physical 
picture, at least, for r and z directions. 

Shear viscosity: Let us evaluate the shear viscosity 
from the data of the stress tensor. The theoretical shear 



viscosity for frictionless granular fluids, which depends 
on temperature and volume fraction, is given by 
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with strain rate D rz = (dv r /dz + dv z /dr)/2, rj*(cf>) = 
r] k *[l + 40x(l + e)/5] + 3 7 */5, 7? fe * = [1 - 2(1 + e)(l - 
3e)<f> X /5}/(u;-C/2), 7 * = 1280 2 X (1 + e)(l - c*/32)/57r, 
v* = X [l - (1 - e) 2 /4][l - c*/64], C = 5 X (1 - e 2 )(l + 
3c732)/12and c* = 32(l-e)(l-2e 2 )/[81-17e+30e 2 (l- 
e)] M- 
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FIG. 5: (Color online) Scaled VDF f(c a ) for <j/ = 0.90. 
Empty points and filled points denote the frictionless case 
and the frictional case, respectively. The black solid line is 

/>«) ~~ 
exp(- 



= exp(— c a /2) /\/2n and the dashed line is f(c a ) 
V2\c a \)/V2. 



For the frictional case, yield stress cry, which is the 
residual stress without deformation may exist in general. 
Thus, the constitutive equation in Eq. (6) is replaced 
by o~rz = o~y — fjDrz, for this case. However, in this let- 
ter, we assume uy = 0. The reason for the absence of 
the yield stress is based on the following three observa- 
tions. The first reason is the velocity distribution func- 
tion. It is known that the significant effect of Coulombic 
slip may appear in the difference of velocity distribution 
functions (VDF), which is characterized by flatness of 
them. The VDF are near to the Gaussian for friction- 
less cases. On the other hand, the VDF for frictional 
cases are exponential-like [3r| . We scale VDF f(v a ) as 
f(v a ) = v~of(c a ) with / dc a f(c a ) = J dc a c 2 a f{c a ) = 1, 
/ dc a c a f(c a ) = and v a o = \j2T a /m for each a. 
Scaled VDF for each velocity components are shown in 
Fig. [5j where all of the VDF are near to Gaussian 
/(c Q ) = cx P(^ c a/2) /v27T even for the frictional case, be- 
cause friction constant /i = 0.2 may be small, and are far 
from exponential-like VDF f(c a ) = exp(— \f2\c a \) / 
The flatness, which is defined as (x 4 )/(x 2 ) 2 — (x A ) for 
(x 2 ) = 1 with (•••) = Jdxf(x)---, is summarized in 
TABLE I for 4> = 0.90. It should be noted that the flat- 
ness with Gaussian VDF is 3.0 and that with exponential 
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TABLE I: Flatness for the frictionless and the frictional case. 





z 


r 


e 


frictionless case 


2.87 


2.86 


3.41 


frictional case 


2.70 


2.98 


3.71 



TABLE II: Extrapolated yield stress -a Y X 10 3 



^^V^tar 

4>o 


0.30 


0.50 


0.70 


0.90 


0.30 


6.41 ± 3.3 


0.701 ± 1.5 


0.434 ± 0.66 


0.447 ± 0.26 


0.40 


7.54 ± 5.8 


-0.159 ± 3.2 


0.538 ± 1.1 


0.364 ± 0.52 


0.45 


8.79 ± 2.3 


0.800 ± 2.1 


0.440 ± 0.94 


0.783 ± 0.44 


0.50 


6.79 ± 2.8 


-0.574 ± 1.9 


0.632 ± 0.94 


0.608 ± 0.32 


0.80 


7.00 ± 4.8 


2.19 ± 1.7 


-0.953 ± 1.3 


0.392 ± 0.50 


0.90 


5.21 ± 2.7 


-0.540 ± 2.9 


-0.396 ± 1.2 


0.243 ± 0.83 



VDF is 6.0. Although the flatness with VDF for 9 devi- 
ates slightly from 3.0, it is still far from 6.0, and thus, the 
effect of Coulombic slip with friction constant /i = 0.2 is 
not significant. The second reason is the small Coulombic 
constant. In this case, renormalization of restitution co- 
efficient is known to be valid {32H341] . We stress here that 
the residual stress for the frictional granular fluid with 
small Coulombic constant does not exist in the usual 
setup[13, HH even for the jamming transition. More- 
over, once we assume that shear viscosity corresponds 
to the value from the kinetic theory rj(r) — r/ki n (4>, T g ), 
the extrapolated ay are obtained at each mesh, through 
a rz (r) = u Y (r) - rfkin((f>(r),Tg(r))D rz (r). As a result, 
ay — holds for r/i? tar > 0.40 (TABLE II ). Thus, 
ay = is a self consistent assumption if we adopt the 
kinetic theory. For these reasons, we assume cry = 0. 
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FIG. 6: (Color online) Non-dimensional shear viscosity n* , 
which is defined in Eq. (7), for several 4>o/4>{ cc in 0.2i? t ar < 
r < -Rtar- Black solid lines denote theoretical curves. 

Now, let us try to compare the theoretical expres- 
sion in Eq. (7) with the numerical measured shear vis- 
cosity. We estimate strain rate as dv r (r, Az/2)/dz ~ 
(v r (r, 3Az/4) - v r (r, Az/4))/(Az/2) and dv z (r, z)/dr ~ 
(v z (r + Ar/2 7 z) — v z (r — Ar/2, z))/Ar. Since we evaluate 
the physical quantities near the wall, the mesh < z < 



Az is divided into < z < Az/2 and Az/2 < z < Az 
to calculate dv r (r, Az/2)/dz and < r < R tal is divided 
into < r < Ar/2, Ar/2 < r < 3Ar/2, • • • . The compar- 
ison of rf , which is the non-dimensional shear viscosity 
in Eq. (7), for 0.2 < r/i? tar < 1.0 is shown in Fig. [BJ Al- 
though there is a slight deviation between them for large 
(f> i.e. small r, which may be the effect of the singular- 
ity in the center r = 0, the theoretical curve reproduces 
other numerical results. We can, thus, conclude that the 
flow has the finite shear viscosity which has the same 
order of the predicted value by the kinetic theory. The 
reason why the simulations are approximately described 
by the Euler equation is that the strain rate itself is small 
i.e. D rz d/uQ < 4.5 x 10 -2 , and thus a rz /a zz is small. 

Dicussion: There exists large normal stress differences, 
which may not be considered in the kinetic theory. We 
found that the normal stress vertical to the wall is almost 
twice as large as the other diagonal components of the 
stress tensor. A unidirectional flow can be distributed 
into two directions in the usual time revolution. This 
mechanism might be easily understood by the time re- 
versal flow in which two directional flow merge into a 
unidirectional flow. This picture is possible to be consid- 
ered because the dissipation is not crucially important, 
at least, for the scattering angle 27]. Thus, the fluctu- 
ations in r and 9 components may be half of that in z 
component. Therefore, we may understand the relation 
T z ~ 2T g ~ 2T r . 

Conclusion: We have numerically investigated the 
granular jet which impacts on a fixed wall. We revealed 
that the granular flow after the impact has the finite 
shear viscosity which has the same order of the predicted 
value by the kinetic theory, and thus the similarity be- 
tween the granular flow and the perfect fluid is superfi- 
cial, which comes from a small strain rate. This result 
is consistent with the result of the experiment [l(| and 
is contrast to the two dimensional simulations for an in- 
viscid fluid 0, [2l[. We have assumed ay — 0, judging 
from VDF, small Coulombic constant and extrapolated 
cry. This assumption is strong one for the comparison 
between the kinetic theory and the our data, but the va- 
lidity of this assumption for larger Coulombic constant 
will be discussed later. Although both the pressure and 
the viscosity are not far from the predictions by the ki- 
netic theory, there exists a large normal stress difference 
in contrast to the case of the kinetic theory. Our results 
may shed the light on the internal fluid structure under 
a strong nonequilibrium situation, i.e. the impact pro- 
cesses of the granular jet. 
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